Minimal Model for Sand Dunes 
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We propose a minimal model for aeolian sand dunes. It combines an analytical description of 
the turbulent wind velocity field above the dune with a continuum saltation model that allows for 
saturation transients in the sand flux. The model provides a qualitative understanding of important 
features of real dunes, such as their longitudinal shape and aspect ratio, the formation of a slip face, 
the breaking of scale invariance, and the existence of a minimum dune size. 



Sand dunes develop wherever sand is exposed to an 
agitating medium (air, water . . . ) that lifts grains from 
the ground and entrains them into a surface flow. The 
diverse conditions of wind and of sand supply in different 
regions on Earth give rise to a large variety of shapes of 
aeolian dunes |l]-[|. Moreover, dunes have been found 
on the sea-bottom and even on Mars. Despite the long 
history of the subject, the underlying physical mecha- 
nisms of dune formation are still not very well under- 
stood. How are aerodynamics (hydrodynamics) and the 
particular properties of granular matter acting together 
to create dunes? How is the shape of a dune maintained 
when it moves? In the following we propose a "minimal 
model" for aeolian sand dunes to address such questions. 
Although it refers only to rather generic properties of the 
wind velocity field and the laws of aeolian sand transport, 
it can make interesting qualitative predictions that are 
not sensitive to the simplifying assumptions, e.g. about 
the surface profile, the development and position of the 
slip face, dune migration etc. Using results from turbu- 
lent boundary layer calculations J4|-|6|] , we will propose an 
approximate analytical description of the surface shear 
stress exerted by the wind onto a heap of sand. This 
will be combined with a saltation model |?J that allows 
for saturation transients of the sand flux, which are an 
essential element of a consistent description of dunes. 

We start by summarizing some basic knowledge about 
aeolian sand transport and saltation. The mean turbu- 
lent wind velocity above a plane surface increases loga- 
rithmically with height. Its magnitude is specified by a 
characteristic velocity called the "shear velocity" it* and 
defined by ul = tq/p with To the average surface shear 
stress (far away from any obstacle) and p the density of 
air. On a surface covered with sand, the wind entrains 
some grains into a surface layer flow if the shear velocity 
exceeds a threshold value. The grains advance mainly by 
an irregular hopping process, thereby reducing the wind 
velocity in the surface layer. A unique relation between 
the shear stress r and the sand flux q is thus established 
in the equilibrium state. If r is not too close to the 



threshold, one has approximately 
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The index s emphasizes that this simple relation is re- 
stricted to situations where the flux is saturated. Ac- 
cording to Eq. ([!]), the changing wind shear stress above 
a heap of sand h(x,y) is responsible for flux gradients, 
which cause erosion and deposition. Due to mass conser- 
vation, the flux gradient dq/dx on a slice h(x) parallel to 
the wind direction gives rise to a change in height with 
time dh/dt, and thus to migration of the surface profile 
in the wind direction with a velocity 



v cx dq/dh 



(2) 



This velocity is much smaller than the velocity of the 
saltating grains, so that the terrain can be assumed to 
be stationary for considerations concerning the wind and 
saltation dynamics. To obtain q s and the corresponding 
v s = v(q s ) from Eqs.(0) and (||), we thus need to know the 
stationary surface shear stress perturbation f = t/tq — 1 
caused by the heap. This depends on its shape. 

In a first attempt to model f, one might try the affinc 
relation f oc h and combine it with the saturated flux 
approximation q = q s . If this was true, the surface veloc- 
ity v would increase with height (dv/dh > 0) due to the 
nonlinearity of Eq. ([!]). This in turn would decrease the 
upwind slope and increase the lee slope, and thus even- 
tually create a slip face on the lee side, where sand slides 
down in avalanches. Since there is so far nothing in the 
model to stop the further decrease of the upwind slope, 
our model "dune" would then start to decrease in height 
and finally become flat, whatever its windward profile. 
This "zeroth order model" is obviously missing an es- 
sential ingredient. In the following, we will successively 
demonstrate how the modeling of both the shear stress 
and the flux have to be refined to achieve a consistent de- 
scription of dune formation and migration. In particular, 
we will show that a subtle balance of two quantitatively 
small but qualitatively crucial corrections to the zeroth 
order model determine the shape of aeolian heaps and 
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FIG. 1. Surface shear stress t(x)/tq (upper curves) above 
the profiles of Eq.(|^) (lower curves: Jl dashed, /g solid, / c 
dotted line) according to Eq.(§) with aH/L = 0.45, /3 = 0.25. 
The small windward shift of the maximum shear stress with 
respect to the apex is crucial for dune formation. 

dunes. 

An improved description for the shear stress can be ob- 
tained from turbulent boundary layer calculations. The 
surface shear stress perturbation T x _m_ wind direction 
caused by a smooth bump h{x, y) is 



3~ cy\^~x"\ — 
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Fx V {h(x,y)} . 



(3) 



Here, we have neglected a minor term oc In \k x L\, and ab- 
breviated the Fourier transformation from x, y to k x , k y 
by T xy . The parameters a and (3 depend logarithmically 
on the ratio L/zq, where L is the characteristic length of 
h(x, 0) and z is a measure of the surface roughness [ [To] . 
For the following discussion we evaluate Eq. (g) on a slice 
h(x) = h(x, 0) of a transverse dune (/i independent of y). 
We obtain (h' := dh/dx) 



t(x) 
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Minor corrections for dunes with a finite transverse width 
can be calculated jll| but are omitted here for simplicity. 
For certain heap profiles h(x) = Hf(x/L) (Fig. |l|) with 



h = (1 



Ig = eT 



fS = cos" 



(5) 



Eq.([|) can be evaluated analytically. (The cosine 
profile /q is understood to vanish outside the region 
[-7I-/2, 7r/2].) The plots of Eq.(@) in Fig. @ share several 
crucial properties not present in the affine approximation 
f oc h. First, we observe from a dimensional analysis of 
Eq.(||) that f oc H/L for a given shape /. Hence, as one 
expects from the scale invariance of turbulence, the am- 
plification of the shear velocity at the top of the profile is 
determined by its aspect ratio H/L and is essentially in- 
dependent of the absolute height H |l^| . Secondly, at the 
tails of the profiles in Fig. the shear stress decreases 
below its asymptotic value To. This effect is particularly 
pronounced for the profile that has a discontinuous 
second derivative. Further, as a consequence of the sec- 
ond term in Eq. (^J) , the surface shear stress is asymmetric 
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FIG. 2. The surface velocity v s (x) of the cosine-shaped 



heap fc(x) (gray) according to Eqs.( 



for q = 3.2, 



/3 = 0.25, and varying aspect ratios H/L = 0.01 (dashed), 
0.14 (solid), 0.27 (dotted). Note the tendency to increase the 
(unstable) lee slope and to drive the windward slope towards 
a (stable) optimum value that depends on a and (3. If fc is 
interpreted as the envelope of a dune and its separation bub- 
ble, the slip face must be located near the sharp drop of v s {x) 
slightly upwind from the top of the envelope. 



even for symmetric profiles like those of Eq.(|J). For 
q = q s , there can thus be a net deposition on a symmetric 
heap of sand and the shift of the shear stress maximum 
with respect to the top of the heap [|l2| allows deposition 
at the top of the heap. An initially flat heap of sand can 
spontaneously grow and increase its windward slope. 

To understand this instability of a plane sand surface, 
the migration velocity v s (x) of a cosine shaped heap of 
sand fc(x) is depicted in Fig. ||. The decrease of v s (x) on 
the lee side reveals the anticipated self-amplifying insta- 
bility of the lee slope, which gives rise to the formation 
of the slip face. More interestingly, Eq.(^) renders v s (x) 
approximately constant over almost the whole windward 
side if H/L is close to a stable optimum value. Slightly 
better results can be achieved for slightly lower n (with 
slightly larger H/L) but not for the profiles fc and Jl, for 
which v s is always non-uniform. The dashed and dotted 
lines in Fig. || were obtained for a smaller and a larger 
aspect ratio and represent a steepening (v' s < 0) and 
flattening (v' s > 0) of the windward side, respectively. 
Altogether, Fig. || suggests that the coupled Eqs.(|l|) and 
(||) drive a heap of sand towards a "dune" with a cosine- 
like windward profile of a preferred aspect ratio, and a 
slip face on the lee side. 

From the preceding qualitative analysis it is not yet 
obvious that this process converges to a translation in- 
variant steady state solution. Several previous studies 
using similar descriptions either did not scrutinize the 
long time behavior of their models or failed to 

obtain stable dunes Jl5|,[l6). Indeed, the present model 
is still insufficient. The problem with the quasi-linear 
Eqs.(||) and (Q) is that they were derived for smooth hills 
and cannot account for flow separation. Once a slip face 
has developed, with a slope of about 32° — 35° that ter- 
minates in a sharp brink, they can no longer be applied. 
As in the textbook example of a backward facing step, 
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flow separation occurs at the brink and large eddies de- 
velop in the wake region. The recirculating flow in this 
"separation bubble" is bounded by a separating stream- 
line against the region of laminar flow. Following a sug- 
gestion by Zeman and Jensen |l5| , we assume that in 
order to calculate the shear stress on the windward side, 
Eqs.(|) and (|) may be applied to the envelope of a dune 
and its separation bubble (inset of Fig. [|) in place of the 
actual surface profile of the dune. In our numerical com- 
putations we model this effective separation bubble by a 
third order polynomial that is the smooth continuation 
of the windward dune profile with a maximum slope of 

14° g. 

If this generalization is accepted, and if — for qualita- 
tive purposes — one identifies the profile with such 
an envelope, one can draw further interesting conclusions 
from Fig. [2| First, the constant v s (x) obtained for the 
windward side of the cosineprofiles i/g (n ~ 2), provides 
strong evidence that Eqs.(P and (||) allow a dune with 
a cosine-like windward side of suitable aspect ratio to 
remain invariant upon translation by aeolian sand trans- 
port. The steep decrease of v s (x) at a certain distance 
upwind from the top of the envelope, which is a conse- 
quence of the maximum of t(x) being shifted with respect 
to the top of the envelope, must then coincide with the 
position of the brink. Secondly, the stable dune shape 
is scale invariant |l(J . Together with the scaling relation 
f cx H/L mentioned above, this implies that the migra- 
tion velocity v of dunes obeys v cx v?JL, in accord with 
the empirical observation of a roughly reciprocal relation 
between dune velocity and dune height 

However, one has to bear in mind that Eq.(|l|) is re- 
stricted to a completely saturated saltation layer. Sat- 
uration transients in the sand flux, as they occur under 
variable wind or sand conditions, introduce a new charac- 
teristic length, the saturation length £ s jjj related to (but 
distinct from) the mean saltation length of the grains. 
Shape invariance is therefore restricted to large dunes 
(L/£ s — > oo). Estimates [Q of typical absolute values of 
£ s suggest that real dunes are not necessarily large in this 
sense, in accord with field measurements |18j that report 
size dependent dune shapes. Eq.(|l|) also fails at bound- 
aries sand/ground. It predicts deposition at the wind- 
ward foot of an isolated dune, where the shear stress de- 
creases (Fig. ||). Previous numerical studies tried 
to avoid this by focusing onto the short time behavior 
and by ad hoc heuristic methods such as a "smoothing 
operator" jl3| or an "adaptation length" [ fbif . 

In the following, we demonstrate how these problems 
are naturally cured if Eq. ([!]) is replaced by a slightly more 
elaborate description derived from a continuum saltation 
model 0. In the steady state (d/dt ~ 0), that model 
reduces to the differential equation 




£ s dq/dx = q(l- q/q s ) 



(G) 



for the sand flux q, and explicit expressions for £ s (t) and 
q s (r). For wind velocities well above the entrainment 
threshold, £ s (t) levels off, and the saturated flux q s 



FIG. 3. Steady state dune shapes for various dune lengths 
(roughly 50 ^ s . . . 100 obtained from Gaussian heaps by 
integrating Eqs.(|2|) and (^). Flow separation at the brink 
was modeled heuristically by introducing a separation bub- 
ble (dashed lines in the inset). Periodic boundary conditions 
were applied to conserve the sand mass and thus allow for a 
steady state. The influx/outflux vanishes for the dunes with 
slip face but is finite for the heap. The time evolution of the 
latter is shown in Fig 0b. (Plot rescaled for visualization.) 



converges to Eq. ([!]). The detailed forms used in our 
numerical integrations, which involve the entrainment 
threshold, can be found in Q but are not important here. 
During the integration of Eqs.(|J) and (Q), the lee slope 
of the heap is constantly evaluated. If it exceeds 14°, 
the separation bubble is introduced for the shear stress 
calculation as described above. Additionally, the surface 
shear stress is set to zero within the separation zone, and 
surface avalanches are induced if the slope exceeds the an- 
gle of repose (34°). For defimteness we choose the values 
a = 3.2, f3 = 0.25 for the parameters of Eq.([|). Start- 
ing from a Gaussian heap, we obtain the steady state 
profiles shown in Fig. [| Their general shape and av- 
erage windward slope of about 10° are in good accord 
with recent field measurements |l8|, 20| . The shape of 
the largest dune (indeed well described by /g, n « 2) 
nicely confirms our expectations from Fig. |^, whereas for 
small heaps, the instability of the lee slope is washed out 
by the saturation transients. The upwind shift of t(x) 
with respect to h(x) (which is of the order of (3L) and 
the lag of q(x) with respect to t(x) and q s (x) (of the 
order of £ s ) compete to establish a critical heap size for 
slip face formation. Smaller heaps do not spontaneously 
develop a slip face, while somewhat larger heaps evolve 
into a steady state dune with a windward slope that is 
lowered and a slip face that is shifted downwind com- 
pared to the asymptotic shape. The critical dune size 
depends on the wind velocity through (and scales as) the 
saturation length £ s (t) « £ s (tq) (« for large r). 

Fig. [| shows the growth histories of two heaps of dif- 
ferent mass and initial aspect ratio H/L. It demonstrates 
the spontaneous formation of a slip face that is stable for 
the larger heap but finally turns out to be unstable for 
the smaller one. For intermediate sizes, flow separation 
can stabilize an existing slip face although it would not 
have developed from a flat configuration. This is the case 
for the two smallest dunes with slip face in Fig. ||, which 
were produced from a steep initial configuration (as in 
Fig. |^b) . As well as the critical heap size for slip face 
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FIG. 4. Growth histories for two Gaussian heaps of differ- 
ent mass and initial aspect ratio (periodic boundary condi- 
tions). The steep initial configuration of the heap in part b) 
gives rise to the temporary formation of a slip face, which 
turns out to be unstable. Note that the distances to reach 
the steady state shape are different. (Plots rescaled for visu- 
alization.) 



formation (out of a flat initial condition), there exists 
thus a somewhat smaller critical dune size for slip face 
destruction and an intermediate regime, where the pres- 
ence of a slip face depends on the initial conditions. This 
implies, in particular, that a slip face that develops out of 
a flat initial condition has always a finite height. Indeed, 
we find that slip face destruction also occurs at a finite 
height. 

In summary, we proposed a minimal model for aeolian 
sand dunes. The advantage of the simple Eqs. (p[) and 
(|J) over more elaborate nonlinear models 0-|23] for the 
wind shear stress is that they combine a high compu- 
tational efficiency with good accuracy ||, and can thus 
be applied to simulate large scale desert topographies. 
The drawback is that we had to introduce the heuristic 
separation bubble [jl7||. We pointed out that the satu- 
ration transients contained in Eq.(|^) break the expected 
scale invariance and are responsible for the existence of 
a critical heap size for slip face formation and a critical 
dune size for slip face destruction. Further, we analyzed 
the mechanisms responsible for the longitudinal shape 
and aspect ratio of dunes and heaps, and for slip face 
formation. It is an interesting question whether (and 
how) these results can be generalized to full 3d prob- 
lems. The deviations of the wind velocity from the main 
wind direction are small above an obstacle that is flat in 
the transverse direction. One can thus expect that the 
present approach applies to transverse dunes and, with 
minor corrections to Eq.(Q) pl|j, in the vicinity of the 
symmetry plane of a crescent-shaped barchan dune. As 
a direct consequence of our above results, only 3d heaps 
larger than a certain minimum size can spontaneously 
develop a slip face, which then has a finite height and 
width. This minimum size should scale as £ s , which is 
a decreasing function of w» H], Preliminary results from 
an extension of the model to 3d support this picture. 
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